Accompanying readings: Chapters 4, 5
In [2]:
%pylab inline
In [3]:
import nengo
import numpy
n = nengo.neurons.LIFRate()
x = numpy.zeros(100)
for i in range(10):
x[i*10:(i+1)*10]=numpy.random.uniform(-1,1)
figure()
title('$x$ value changing over time')
plot(x)
ylabel('$x$')
ylim(-1,1)
xlabel('Time (ms)')
figure()
title('Firing rate changing over time')
plot(n.rates(x, gain=1.5, bias=1))
ylabel('Firing rate (Hz)')
xlabel('Time (ms)')
show()
Bilipid cell membrane acts as a capacitor
Flow of ions through the cell membrane ion channels is a resistor
Typical values
For lots more info, see Unit 2 in this online course
Kirchhoff's Current Law
Ohm's Law
Capacitor
Applying these laws to RC circuits:
$ \begin{align} V(t) = \frac{R}{\tau_{RC}} \int_0^t e^{-(t-t')/\tau_{RC}} J_M(t')\; dt' \end{align} $
I.e. the voltage right now (at $t$) depends on all past current input, $J_M(t')$, where each input is weighted by a function that exponentially decays as it gets further away (in the past) from the current time (a kind of memory).
If $J_M$ is constant:
$ \begin{align} V(t) &= \frac{R J_M}{\tau^{RC}} e^{-t/\tau_{RC}}\int_0^t e^{t'/\tau_{RC}} \; dt' \\ &= J_M R (1-e^{-t/\tau_{RC}}) \end{align} $
How long does it take to get to threshold?
How long between spikes?
Firing rate
Simplify by assuming $V_{th}=1$ and $R=1$
But we can't assume a constant $J_M$ in a real brain!
In [71]:
import nengo
from nengo.utils.matplotlib import rasterplot
model = nengo.Network(label='Spiking Neurons')
with model:
stimulus = nengo.Node(0)
ens = nengo.Ensemble(1, dimensions=1,
encoders = [[1]],
intercepts = [-.5],
max_rates= [100])
nengo.Connection(stimulus, ens)
spikes_p = nengo.Probe(ens.neurons, 'spikes')
voltage_p = nengo.Probe(ens.neurons, 'voltage')
stim_p = nengo.Probe(stimulus)
sim = nengo.Simulator(model)
sim.run(.5)
t = sim.trange()
figure(figsize=(6, 3))
ax = gca()
ax.plot(t, sim.data[voltage_p],'g')
rasterplot(t, sim.data[spikes_p], ax=ax)
ylim((0,2))
ylabel("Voltage")
xlabel("Time");
In [72]:
import numpy as np
from nengo.utils.functions import piecewise
x = dict(zip(np.linspace(0,.1,10), np.random.uniform(-1,1,10)))
with model:
stimulus.output = piecewise(x)
sim = nengo.Simulator(model)
sim.run(.1)
t = sim.trange()
figure()
title('$x$ value changing over time')
plot(t, sim.data[stim_p])
ylabel('$x$')
ylim(-1,1)
xlabel('Time (s)')
figure()
title('Firing rate changing over time')
plot(ens.neuron_type.rates(sim.data[stim_p], gain=sim.data[ens].gain, bias=sim.data[ens].bias))
ylabel('Firing rate (Hz)')
xlabel('Time (ms)')
figure(figsize=(6, 3))
ax = gca()
rasterplot(t, sim.data[spikes_p], ax=ax)
ax.plot(t, sim.data[voltage_p],'g') #, drawstyle='steps-post') #if want to see exact timing
ylim((0,2))
title('Voltage and spikes over time')
ylabel("Voltage")
xlabel("Time (s)");
In [80]:
#have to run previous cells
with model:
stimulus.output = lambda t: np.sin(10*t)
sim = nengo.Simulator(model)
sim.run(.6)
nengo.utils.matplotlib.plot_tuning_curves(ens, sim)
t = sim.trange()
figure(figsize=(10,4))
plot(t, sim.data[stim_p], label='x')
ax = gca()
ax.plot(t, sim.data[voltage_p],'g', label='v')
ylim((-1,2))
ylabel('$x$')
xlabel('time (s)')
legend(loc='lower left');
rasterplot(t, sim.data[spikes_p], ax=ax.twinx(), use_eventplot=True)
title('Voltage and spikes over time')
ylabel("Voltage");
xlabel("Time");
ylabel('$neurons$')
ylim((-1,2))
Out[80]:
So that's what happens when we allow $x(t)$ to vary over time
How do decode this?
Can we stick with linear decoding?
In [81]:
#assumes has done %pylab inline
import numpy as np
import nengo
from nengo.utils.matplotlib import rasterplot
model = nengo.Network(label='Two Neurons')
with model:
stim = nengo.Node(lambda t: np.sin(10*t))
ens = nengo.Ensemble(2, dimensions=1,
encoders = [[1],[-1]],
intercepts = [-.5, -.5],
max_rates= [100, 100])
nengo.Connection(stim, ens)
stim_p = nengo.Probe(stim)
spikes_p = nengo.Probe(ens.neurons, 'spikes')
sim = nengo.Simulator(model)
sim.run(.6)
nengo.utils.matplotlib.plot_tuning_curves(ens, sim)
t = sim.trange()
figure(figsize=(12, 6))
ax = gca()
plot(t, sim.data[stim_p],'g')
ylabel("Output")
xlabel("Time");
rasterplot(t, sim.data[spikes_p], ax=ax.twinx(), use_eventplot=True)
axis('tight')
ylabel("Neuron");
In [82]:
import numpy as np
import nengo
model = nengo.Network(label='Decoding Neurons')
with model:
stim = nengo.Node(lambda t: np.sin(10*t))
ens = nengo.Ensemble(2, dimensions=1,
encoders = [[1],[-1]],
intercepts = [-.5, -.5],
max_rates = [100, 100])
temp = nengo.Ensemble(10, dimensions=1)
nengo.Connection(stim, ens)
connection = nengo.Connection(ens, temp) #This is just to generate the decoders
stim_p = nengo.Probe(stim)
spikes_p = nengo.Probe(ens.neurons, 'spikes')
sim = nengo.Simulator(model)
sim.run(.6)
t = sim.trange()
x = sim.data[stim_p][:,0]
A = sim.data[spikes_p]
gamma=np.dot(A.T,A)
upsilon=np.dot(A.T,x)
d = np.dot(np.linalg.pinv(gamma),upsilon)
xhat = np.dot(A, d)
figure(figsize=(8,4))
plot(t, x, label='x')
plot(t, xhat)
ylabel('$x$')
#ylim(-1,1)
xlabel('Time');
In [85]:
import numpy as np
import nengo
from nengo.utils.matplotlib import rasterplot
from nengo.dists import Uniform
model = nengo.Network(label='Decoding Neurons')
N = 50
with model:
stim = nengo.Node(lambda t: np.sin(10*t))
ens = nengo.Ensemble(N, dimensions=1,
max_rates=Uniform(100,200))
temp = nengo.Ensemble(10, dimensions=1)
nengo.Connection(stim, ens)
connection = nengo.Connection(ens, temp) #This is just to generate the decoders
stim_p = nengo.Probe(stim)
spikes_p = nengo.Probe(ens.neurons, 'spikes')
sim = nengo.Simulator(model)
sim.run(.6)
x = sim.data[stim_p][:,0]
A = sim.data[spikes_p]
gamma=np.dot(A.T,A)
upsilon=np.dot(A.T,x)
d = np.dot(np.linalg.pinv(gamma),upsilon)
xhat = np.dot(A, d)
t = sim.trange()
figure(figsize=(12, 6))
ax = gca()
plot(t, sim.data[stim_p],'b')
plot(t, xhat,'g')
ylabel("Output")
xlabel("Time");
rasterplot(t, sim.data[spikes_p], ax=ax.twinx(), use_eventplot=True, color='k')
axis('tight')
ylabel("Neuron");
In [3]:
import numpy as np
dt = 0.001
sigma = 0.007
t_h = np.arange(200)*dt-0.1
h = np.exp(-t_h**2/(2*sigma**2))
h = h/norm(h,1)
figure()
plot(t_h, h)
xlabel('t')
ylabel('h(t)');
In [4]:
import nengo
from nengo.utils.matplotlib import rasterplot
model = nengo.Network(label='Decoding Neurons')
with model:
stim = nengo.Node(lambda t: np.sin(10*t))
ens = nengo.Ensemble(2, dimensions=1,
encoders = [[1],[-1]],
intercepts = [-.5, -.5],
max_rates = [100, 100])
temp = nengo.Ensemble(10, dimensions=1)
nengo.Connection(stim, ens)
connection = nengo.Connection(ens, temp) #This is just to generate the decoders
stim_p = nengo.Probe(stim)
spikes_p = nengo.Probe(ens.neurons, 'spikes')
sim = nengo.Simulator(model)
sim.run(.6)
t = sim.trange()
x = sim.data[stim_p]
d = sim.data[connection].weights.T
fspikes1 = np.convolve(sim.data[spikes_p][:,0], h, mode='same')
fspikes2 = np.convolve(sim.data[spikes_p][:,1], h, mode='same')
A = np.array([fspikes1, fspikes2]).T
xhat = np.dot(A, d)
figure(figsize=(8, 4))
ax = gca()
plot(t, sim.data[stim_p],'g')
ylabel("Output")
xlabel("Time");
rasterplot(t, sim.data[spikes_p], ax=ax.twinx(), use_eventplot=True, color='k')
ylabel("Neuron")
figure(figsize=(8,4))
plot(t, x, label='x')
plot(t, fspikes1*d[0])
plot(t, fspikes2*d[1])
ylabel('$x$')
xlabel('time (s)')
figure(figsize=(8,4))
plot(t, x, label='x')
plot(t, xhat, label='x')
ylabel('$x$')
xlabel('time (s)');
Since there are two neurons and they are opposites of each other, we can assume that $d_1=-d_2$ (i.e. they both contribute equally, but oppositely)
Let's call $a_1(t)-a_2(t)=r(t)$, for response
$\hat{x}(t)=r(t) * h(t)$
In [5]:
import numpy as np
x = np.random.uniform(-1,1,500)
plot(np.arange(500)*0.001, x);
In [6]:
from nengo.processes import WhiteSignal
dt = 0.001
max_freq = 10
model = nengo.Network('White Noise')
with model:
sig = nengo.Node(output=WhiteSignal(1, high=10, rms=0.5))
sig_p = nengo.Probe(sig)
sim = nengo.Simulator(model)
sim.run(1.0)
t = sim.trange()
x = sim.data[sig_p]
figure(figsize=(8, 4))
ax = gca()
plot(t, sim.data[sig_p])
ylabel("Output")
xlabel("Time");
In [11]:
import numpy as np
T=1.0
dt=0.001
signal = sim.data[sig_p][:,0]
sig_freqs = np.fft.fftshift(np.fft.fft(signal))
t = np.arange(int(T/dt))*dt
freq = np.arange(int(T/dt))/T - (T/dt)/2
figure()
plot(freq, np.abs(sig_freqs))
xlim(-500, 500)
xlabel('freq (Hz)')
ylabel('$|X(\omega)|$')
figure()
plot(freq, np.abs(sig_freqs))
xlim(-50, 50)
xlabel('freq (Hz)')
ylabel('$|X(\omega)|$')
figure()
plot(t, signal)
xlabel('time (s)')
ylabel('$x(t)$');
$\hat{X}(\omega) = R(\omega)H(\omega)$
We want to find $H(\omega)$ that minimizes the error
But, we know that our signal $x(t)$ can be written in the frequency domain as a sum of sine waves (that's the Fourier transform)
So, we can write our error in terms of different frequency values $\omega$
$E(\omega) = (X(\omega)-R(\omega)H(\omega))^2$
Now we can take the derivative and set it equal to zero to do the minimization
$H(\omega)= {{X(\omega)R^*(\omega)} \over {|R(\omega)|^2}} $
So now we can find $H(\omega)$ given the Fourier transform of our signal $X(\omega)$ and the Fourier transform of the spiking response $R(\omega)$
In [4]:
import nengo
from nengo.processes import WhiteSignal
dt = 0.001
max_freq = 10
T = 1.0
model = nengo.Network('White Noise')
with model:
in_stim = nengo.Node(output=WhiteSignal(T, high=max_freq, rms=0.5))
ens = nengo.Ensemble(2, dimensions=1,
encoders = [[1],[-1]],
intercepts = [-.3, -.3],
max_rates= [100, 100])
temp = nengo.Ensemble(1, dimensions=1)
nengo.Connection(in_stim, ens)
connection = nengo.Connection(ens, temp) #This is just to generate the decoders
stim_p = nengo.Probe(in_stim)
spikes_p = nengo.Probe(ens.neurons, 'spikes')
dec_p = nengo.Probe(ens)
sim = nengo.Simulator(model)
sim.run(T)
r_plot = (sim.data[spikes_p][:,0]-sim.data[spikes_p][:,1])*dt
r = sim.data[dec_p][:,0] #properly scaled r
signal = sim.data[stim_p][:,0]
t = sim.trange()
figure(figsize=(8, 6))
plot(t, sim.data[stim_p],'g')
plot(t, r_plot,'r')
ylabel("Output")
xlabel("Time")
axis('tight');
In [5]:
freq = np.arange(int(T/dt)) - (T/dt)/2
omega = freq*2*numpy.pi
X = np.fft.fftshift(np.fft.fft(signal))
R = np.fft.fftshift(np.fft.fft(r))
figure()
plot(omega, np.abs(X))
xlabel('$\omega$')
ylabel('$|X(\omega)|$')
xlim(-200,200)
figure()
plot(omega, np.abs(R))
xlabel('$\omega$')
ylabel('$|R(\omega)|$');
Now we can find our optimal filter
$H(\omega)= {{X(\omega)R^*(\omega)} \over {|R(\omega)|^2}} $
In [6]:
Rmag = R*R.conjugate()
Rmag[(Rmag.shape[0]-1)/2 + 1] = (Rmag.shape[0]-1)/2 #set zero DC term
H = X*R.conjugate() / Rmag
h = np.fft.fftshift(np.fft.ifft(np.fft.ifftshift(H))).real
figure()
plot(omega, np.abs(X*R.conjugate()))
xlim(-200,200)
figure()
plot(omega, np.abs(Rmag))
figure()
plot(omega, np.abs(H))
xlabel('$\omega$')
ylabel('$|H(\omega)|$')
xlim(-200,200)
figure()
plot(freq*dt,h)
xlabel('t')
ylabel('h(t)')
xlim(-.2,.2);
In [8]:
fspikes1 = np.convolve(sim.data[spikes_p][:,0], h, mode='same')
fspikes2 = np.convolve(sim.data[spikes_p][:,1], h, mode='same')
A = np.array([fspikes1, fspikes2]).T
d = sim.data[connection].weights.T
xhat = np.dot(A, d)[:,0]
figure()
plot(t, signal)
plot(t, xhat)
show()
print 'RMSE:', np.average((signal-xhat)**2)
In [13]:
sim = nengo.Simulator(model)
sim.run(T)
sig = sim.data[stim_p][:,0]
fspikes1 = np.convolve(sim.data[spikes_p][:,0], h, mode='same')
fspikes2 = np.convolve(sim.data[spikes_p][:,1], h, mode='same')
A = np.array([fspikes1, fspikes2]).T
d = sim.data[connection].weights.T
xhat = np.dot(A, d)[:,0]
t = sim.trange()
figure()
plot(t, sig)
plot(t, xhat)
show()
print 'RMSE:', np.average((sig-xhat)**2)
In [14]:
sig2 = np.zeros_like(sig)
r2 = np.zeros_like(r)
sig2[400:600] = sig[100:300]
r2[400:600] = r[100:300]
figure()
plot(t, sig)
plot(t, r)
vlines(0.1, -10, 10, linewidth=3)
vlines(0.3, -10, 10, linewidth=3)
figure()
plot(t, sig2)
plot(t, r2);
In [15]:
sigma = 0.06 #This should be ~3x longer than the neuron memory
window = np.exp(-(t-0.5)**2/(sigma**2))
sig2w = np.roll(sig, 300)*window
r2w = np.roll(r, 300)*window
figure()
plot(t, sig)
plot(t, r)
vlines(0.1, -10, 10, linewidth=3)
vlines(0.3, -10, 10, linewidth=3)
figure()
plot(t, sig2w)
plot(t, r2w)
show()
In [16]:
# original code ignoring DC
H = (X*R.conjugate()) / (R*R.conjugate())
# new code
sigma_t = 0.025
W2 = np.exp(-omega**2*sigma_t**2)
H = (np.convolve(X*R.conjugate(),W2, 'same')) / (np.convolve(R*R.conjugate(),W2, 'same'))
How well does this work?
In [33]:
in_stim.output.high=10
sim = nengo.Simulator(model)
sim.run(T)
sig1 = sim.data[stim_p][:,0]
r = sim.data[dec_p][:,0] #properly scaled r
X = np.fft.fftshift(np.fft.fft(sig1))
R = np.fft.fftshift(np.fft.fft(r))
sigma_t = 0.06
W2 = np.exp(-omega**2*sigma_t**2)
Rmag = R*R.conjugate()
Rmag[(Rmag.shape[0]-1)/2 + 1] = .1 #remove zero DC term
H = (np.convolve(X*R.conjugate(),W2, 'same')) / (np.convolve(Rmag,W2, 'same'))
h = np.fft.fftshift(np.fft.ifft(np.fft.ifftshift(H))).real
In [34]:
figure()
plot(omega, np.abs(H))
xlabel('$\omega$')
ylabel('$|H(\omega)|$')
xlim(-500, 500)
figure()
plot(t-T/2, h)
xlabel('t')
ylabel('h(t)')
figure()
plot(t-T/2, h)
xlabel('t')
ylabel('h(t)')
xlim(-0.2, 0.2)
Out[34]:
In [35]:
spikes1 = sim.data[spikes_p][:,0]
fspikes1_1 = np.convolve(sim.data[spikes_p][:,0], h, mode='same')
fspikes2_1 = np.convolve(sim.data[spikes_p][:,1], h, mode='same')
spikes_1 = (sim.data[spikes_p][:,0]-sim.data[spikes_p][:,1])*sim.dt
sim = nengo.Simulator(model)
sim.run(T)
sig2 = sim.data[stim_p][:,0]
fspikes1_2 = np.convolve(sim.data[spikes_p][:,0], h, mode='same')
fspikes2_2 = np.convolve(sim.data[spikes_p][:,1], h, mode='same')
spikes_2 = (sim.data[spikes_p][:,0]-sim.data[spikes_p][:,1])*sim.dt
A1 = numpy.array([fspikes1_1, fspikes2_1]).T
A2 = numpy.array([fspikes1_2, fspikes2_2]).T
d = sim.data[connection].weights.T
x1hat = np.dot(A1, d)[:,0]
x2hat = np.dot(A2, d)[:,0]
t = sim.trange()
figure()
plot(t, sig1)
plot(t, x1hat)
plot(t, spikes_1, 'k.')
ylim(-1.5,1.5)
show()
print 'RMSE:', np.average((sig1-x1hat)**2)
figure()
plot(t, sig2)
plot(t, x2hat)
plot(t, spikes_2, 'k.')
ylim(-1.5,1.5)
print 'RMSE:', np.average((sig2-x2hat)**2)
In [36]:
in_stim.output.high=2
sim = nengo.Simulator(model)
sim.run(T)
sig3 = sim.data[stim_p][:,0]
fspikes1_3 = np.convolve(sim.data[spikes_p][:,0], h, mode='same')
fspikes2_3 = np.convolve(sim.data[spikes_p][:,1], h, mode='same')
spikes_3 = (sim.data[spikes_p][:,0]-sim.data[spikes_p][:,1])*sim.dt
A3 = np.array([fspikes1_3, fspikes2_3]).T
d = sim.data[connection].weights.T
x3hat = np.dot(A3, d)[:,0]
t = sim.trange()
figure()
plot(t, sig3)
plot(t, x3hat)
plot(t, spikes_3, 'k.')
ylim(-1.5,1.5)
print 'RMSE:', np.average((sig3-x3hat)**2)
In [37]:
figure()
plot(t-T/2, h)
xlabel('t')
ylabel('h(t)')
xlim(-0.2, 0.2)
Out[37]:
In [38]:
figure(figsize=(8,4))
vlines(dt*np.where(spikes1>0)[0], -1, 1)
plot(t, fspikes1_1*d[0], label='filtered spikes')
ylabel('$x$')
ylim(-1,1)
xlim(0.2, 0.4)
xlabel('time (s)')
legend();
(image stolen from here)
This is the standard simple model
For more complex models, some people use
In [39]:
dt = 0.001
tau = 0.05
t_h = np.arange(1000)*dt-0.5
h = np.exp(-t_h/tau)
h[np.where(t_h<0)]=0
h = h/norm(h,1)
figure()
plot(t_h, h)
xlabel('t')
ylabel('h(t)');
In [43]:
in_stim.output.high=10
sim = nengo.Simulator(model)
sim.run(T)
sig = sim.data[stim_p][:,0]
fspikes1 = np.convolve(sim.data[spikes_p][:,0], h, mode='same')
fspikes2 = np.convolve(sim.data[spikes_p][:,1], h, mode='same')
A = np.array([fspikes1, fspikes2]).T
d = sim.data[connection].weights.T
xhat = np.dot(A, d)
t = sim.trange()
from nengo.utils.matplotlib import rasterplot
figure(figsize=(8, 4))
ax = gca()
plot(t, sig,'g')
ylabel("Output")
xlabel("Time");
rasterplot(t, sim.data[spikes_p], ax=ax.twinx(), use_eventplot=True, color='k')
ylabel("Neuron")
figure(figsize=(8,4))
plot(t, sig, label='x')
plot(t, fspikes1*d[0])
plot(t, fspikes2*d[1])
ylabel('$x$')
ylim(-1,1)
xlabel('time (s)')
figure(figsize=(8,4))
plot(t, sig, label='x')
plot(t, xhat, label='x')
ylabel('$x$')
ylim(-1,1)
xlabel('time (s)');
In [45]:
from nengo.dists import Uniform
with model:
ens = nengo.Ensemble(40, dimensions=1, max_rates=Uniform(100,200))
nengo.Connection(in_stim, ens)
dec_p = nengo.Probe(ens, synapse = 0.005)
sim = nengo.Simulator(model)
sim.run(T)
sig = sim.data[stim_p][:,0]
t = sim.trange()
figure(figsize=(8,4))
plot(t,sim.data[dec_p])
plot(t,sig);
It acts like an infinite amount of training with constant inputs
The filtered spiking neuron model gives some sort of variability around the ideal $a$ value